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ABSTRACT 

In this paper we carry out a quantitative analysis of the three-body systems and map them 
as a function of decaying time and initial configuration, look at this problem as an example 
of a simple deterministic system and ask to what extent the orbits are really predictable. 
We have investigated the behaviour of about 200 000 general Newtonian three-body systems 
using the simplest initial conditions. Within our resolution these cover all the possible states 
where the objects are initially at rest and have no angular momentum. We have determined 
the decay time-scales of the triple systems and show that the distribution of this parameter is 
fractal in appearance. Some areas that appear stable on large scales exhibit very narrow strips 
of instability and the overall pattern, dominated by resonances, reminds us of a traditional 
Maasai warrior shield. Also an attempt is made to recover the original starting configuration 
of the three bodies by backward integration. We find there are instances where the evolution to 
the future and to the past lead to different orbits, in spite of time symmetric initial conditions. 
This implies that even in simple deterministic systems there exists an arrow of time. 
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1 INTRODUCTION 


Dynamical interaction of multiple objects is common in physics 
and astronomy. It occurs over a huge range of scales from atoms in 
molecules to galaxies. If we know the interaction potential, then it 
is possible in principle to calculate the evolution of such a system. 
As is well known, systems of this kind tend to be chaotic, and so is 
the solution of the three-body problem (Prigogine & Stengers 1984; 
Aarseth et al. 1994). 

The two-body problem in a Newtonian gravitational potential is 
simple enough to be solvable in a closed form. Adding a third body 
to the system complicates matters by producing, in general, non- 
analytic solutions. Over a century ago Poincaré (1892) realized the 
chaotic nature of the few body system, but only with the advance of 
computer power and more sophisticated algorithms has it been pos- 
sible to give more quantitative statements in the three-body problem 
(Heggie 1975; Aarseth et al. 1994; Mikkola 1994). 

Here we introduce a simple way to illustrate and study the frac- 
tal nature of the general three-body system. The method we de- 
scribe provides an important tool for visual and qualitative studies of 
the three-body problem. It allows us to handle and analyse all pos- 
sible three-body configurations within our computational accuracy. 
The results visualized in this manner do not suffer of projection ef- 
fects. With this method it is also easy to extend our study to include 
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the influence of different initial conditions such as different veloci- 
ties or masses. These extensions will help us understand the chaotic 
nature of the three-body systems both in much greater detail and in 
a more comprehensive way. We also discuss shortly the numerical 
sensitivity of general three-body problem and show that this may 
indicate the presence of the arrow of time within these systems. 


2 DECAY TIME-SCALE AND FRACTALITY 


We have considered the simplest initial state in an equal mass three- 
body system under Newtonian gravitation. The calculations start at 
rest and have no angular momentum. We have used in our calcula- 
tions a fully chain-regularized Bulirsch—Stoer integrator (Mikkola 
& Aarseth 1990, 1993, 1996). 

To visualize our simulations we use an extended version of a ho- 
mology map (Agekian & Anosova 1967). In this coordinate system 
one body resides initially at (0,0), the second one is at (1,0) and 
the third one is at (x,y). A shield-shaped area contains all possible 
values of (x,y), spanning in abscissa from 0.0 to 1.0 and in ordinate 
from —/3/2 to /3/2 ~ 0.866. A three-body configuration can 
thus be unambiguously specified by a single point in this map. The 
unit time-scale of our (G = 1) system is t = (1/27) (R/1 au)*/? 
(M /Mo)"/ 2 yr, where R is the longest separation of the initial 
triangle and M is the mass of the objects. 

Except for a few double or triple collision points (Umehara & 
Tanikawa 2000), all our three-body systems will eventually breakup. 
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The escaping body will have a non-negative total orbital energy. To 
define the breakup time of the system we also require that the 
escaping body must be at a distance of r > 1 from both remaining 
bodies. This is necessary because during strong interactions the total 
energy of a body may, as expected, become momentarily highly 
positive. 

Because the escaping body often has only a mildly positive 
energy and because the quantity we measure is time (the conju- 
gate coordinate of energy in phase space), our measurement is in 
some respect analogous to a phase-space cross-section of a triple 
system. 

The lifetime of the system may now be plotted as a function of 
the initial configuration. A complicated structure emerges (Fig. 1). 


O 0:3 


Figure 1. The lifetime of gravitational three-body systems containing all 
possible of orbits that were initially at rest. The initial configuration of the 
three-body system is such that one star is located at (0,0), the second one at 
(1,0) and the third one at a coordinate (x,y). The colour pixel at (x,y) shows 
then the time elapsed from the initial configuration to the disintegration of 
the system under Newtonian gravitation. Red denotes short lived systems, 
decaying within one crossing time, orange systems decaying before two 
crossing times etc. The dark blue colours represent systems with the longest 
lifespan (over 13 initial crossing times). The pixel size is 1.25 x 107° . The 
image of the complete shield is symmetric, because we start with all three 
objects at rest. Note the wide resonant bands causing a rapid disintegration 
of the system. Also note how narrower bands abound between wide bands. 


As a curiosity, the overall pattern shows a nice resemblance with 
a traditional African Maasai warrior shield (see e.g. Huntingford 
1961). Since we consider an equal mass system, there is a left-right 
symmetry in the diagram. The non-zero initial velocity, e.g. for the 
third body, will in general break the up-down symmetry. 


2.1 Resonances 


On the largest scales the following regions of instability can be 
identified. A configuration initially close to an isosceles triangle 
(1:1 resonance) is generally highly unstable (the vertical red band 
close to the vertical symmetry line of the map). Systems close to 
odd-digit orbital resonances (1:3, 1:5, 1:7 etc.) are highly unstable, 
because a very close approach takes place at the first close triple 
encounter. These are the bands growing denser towards the left- 
and the right-hand corners of the map. Note also the branching of 
the 1:7 and higher order resonances. Furthermore, resonances 1:2, 
1:4 etc. are completely missing, because the phase in these cases 
is not favourable at the first close encounter to a rapid breakup of 
the system (Saslaw, Valtonen & Aarseth 1974). Systems breaking 
up due to the second close interaction show up as long orange arcs. 
The family of these arcs starts from the central part of the map 
and extends towards the edges of the map. A second system of 
arcs between the previous set and the 1:3 resonance is also due to 
the odd resonances at the second encounter. A band due to similar 
resonances (nearly vertical bands) can be seen at the edge of the 
1:1 resonance. 

All the major resonances appear to avoid a sea of apparent tran- 
quillity around the surroundings of (x, y) = (0.7, 0.5) and the cor- 
responding points in the other quarters of the shield. 


2.2 Local structure 


Locally the boundaries of these resonance regions are sharp from 
the inside. Approaching the boundary from outside one encounters 
ever denser loops and stripes. This structure appears self-similar 
and locally fractal. 

Areas of higher stability are also present in Fig. 1. These fall into 
two categories. The first set are isolated deep trenches on large scales 
at about (0.70, 0.64) and (0.65,0.15) and at respective symmetric 
points of the shield. The second set better visible on small scales 
are configurations just outside some instability regions e.g. around 
the upper branch of the 1:7 resonance at (0.90, 0.23) or below 
the yellow eye at (0.75,0.60) and the corresponding symmetric 
points. 

Figs 2 and 3 show progressive levels of zooming into these areas. 
Fig. 3 represents a 5 x 5 pixel area in the lower left-hand corner of 
Fig. 2. The width of the strips seen in Fig. 3 is comparable to the 
resolution element or 3.125 x 10~°. These stripes are more easily 
seen if the image is viewed at a small angle. It is clear that a tiny 
change in the initial conditions may cause a substantial difference 
in the lifetime and in the evolution of the system. Further details of 
the orbital behaviour within these areas are beyond the scope of this 
paper, except that we may state that these configurations appear to 
be more stable. 


2.3 Global fractality 


To investigate the ratio of stable (locally homogeneous) regions to 
the chaotic regions we compared f,, the ejection time of an orbit, 
tO tmea, the median of eight nearby orbits located on a rectangular 
grid with total width of 3.75 x 1073. The eight configurations are 
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Figure 2. The lifetime of gravitational three-body systems, a zoomed in 
image of Fig. 1. The pixel size here is 2.5 x 1074. Colours as in Fig. 1. 
Note how this region appears quite ‘stable’ in Fig. 1, yet in this image it 
has broken into hundreds of narrow bands of instability. Note that the large 
yellow area breaks into three major parts. This is a universal feature for 
the loops. They correspond to the different object being ejected from the 
system. The point where these three parts meet is a triple collision point. 
The boundary of the yellow area is sharp from inside. If one approaches the 
boundary from the outside one encounters a number of orbits growing in 
density and stability as one approaches the sharp boundary. 


separated by € = +1.25 x 107° in each of the x and y coordinates 
from the initial point. 

If |tmea — tol < 0.5, we considered the point to have a lo- 
cally homogeneous neighbourhood. Otherwise it was considered 
to have high sensitivity to initial conditions, a characteristic of 
chaos. Of the 21 699 cases compared in this manner 6212 or pe = 
28.6 per cent were in locally homogeneous. The remaining 71.4 per 
cent of positions turned out to be highly sensitive to initial condi- 
tions at this resolution implying that on the scale of € = 1.25 x 
107° the fractal area covers about 71.4 per cent of all the systems. 
The implication of this is that if one is able to determine the initial 
state of a triple system to an accuracy of €, one will, with the stated 
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Figure 3. The lifetime of gravitational three-body systems, a zoomed image 
of Fig. 2 showing the smaller details. The pixel size here is 3.125 x 1076. 
The colour scale is linear and is extended to 20 crossing times. The dark 
blue colours show orbits that have not disintegrated in 20 crossing times. 
The light green band in the right upper corner of this image corresponds to 
a disintegration in about seven crossing times. 


probability, be unable to tell the time-scale of the ejection of the 
third body to within half a crossing time. 

If our image is fractal one should find more and more (local) 
order as one decreases the resolution element. This will be difficult 
to do due to restrictions in computation time. We can increase the 
resolution and see if the amount of local homogeneity decreases 
in a self-similar way. For bin sizes one, two, three, four, five and 
10 times the original one we find a decrease in local homogeneity, 
De, 28.6, 16.8, 14.9, 13.8, 13.5 and 8.6 per cent, respectively. To 
within measurement errors pe x €~°°>. We find the same exponent 
if we study the detailed areas shown in Figs 2 and 3. The Hausdorff 
dimension (Hausdorff 1919) of our image is thus Dy = 1.5. 


3 TIME REVERSAL AND THE ARROW 
OF TIME 


The initial conditions we choose for our systems have the important 
mathematical property that the orbits we calculate are time symmet- 
ric with respect to the origin. This implies that the triple systems we 
calculate are transient in the sense that the body expelled must have 
been once captured by the very same bodies forming the final binary 
system (and even in the very same phase but with an opposite sign!). 
Physically this is not quite true, because even a small disturbance to 
the system would cause the system to loose this property. This also 
implies that a simple time reversal at the time of escape will not be 
a useful indicator of the errors involved in these calculations. 

It is possible that the initial conditions are not recovered by two 
reasons: the calculations may not be accurate enough, and errors ac- 
cumulate to the extent where the initial conditions are missed on the 
way back. The second possibility is that the calculation is accurate 
enough, but the individual system is so sensitive that it is practically 
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impossible to find the initial conditions again. See Aarseth et al. 
(1994) for a discussion of predictable and non-predictable orbits. 

To ensure independence from the method our calculations were 
done with two different methods. The method uses a code con- 
sisting of a Bulirsch-Stoer integrator with the chain regulariza- 
tion algorithm (Mikkola & Aarseth 1993). The second code uses a 
simple leapfrog integrator to logarithmic Hamiltonian (Mikkola & 
Tanikawa 1999). Results obtained are in agreement but in analysis 
we prefer the latter code because of its constant step size. 

The orbits were calculated in few stages. First, a forward inte- 
gration was carried out until one of the bodies escapes. Here one 
should notice that criterion for recognizing an escape is not trivial. 
An escape is suspected if the ratio between the semimajor axis of 
one body orbit and the semimajor axis of the binary (i.e. two bodies 
close to each other) is large enough. We have selected this ratio to 
be 100. In addition to this requirement the system has to fulfil one 
of the two following conditions. An escape is identified if the third 
body is on a hyperbolic orbit and is moving away from the binary. 
Alternatively, if one of the bodies gets so far that the ratio exceeds 
>2000 we consider this an escape, even though a proper hyperbolic 
orbit is not present. These escape criteria follow Anosova (1986), 
which includes also a list of different possible escape criteria. We 
have also tried smaller values for the escape criteria with no apparent 
effect on our results. 

At the second stage the escape velocities were reversed and a 
reverse orbit was calculated until the system reached its initial con- 
ditions again. After that the simulation has carried on until the 
system breaks up again. 

The definition of when the initial conditions have been recovered 
depends on the accuracy of the computer and software. We find that 
by using double precision we can expect results to be correct at 
least with accuracy ~10~’. To avoid rounding errors we consider 
only five significant digits. Therefore, we require that the separation 
between original initial state and its recovered counterpart must be 
less than 6 = +5.0 x 1075. The choice of this limit affects naturally 
the ratio of the recovered and unrecovered orbits. The ratio is shown 
in Fig. 4 where different values of 5 have been applied to the data 
set. 

If one would have an infinite accuracy one would obtain an orbit 
where the motion is exactly the same before and after the free-fall 
moment — in other words the orbit is time symmetric. This is a fact 
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Figure 4. The number fraction of recovered and unrecovered orbits as a 
function of required accuracy ô. At around ô = 107% one come to the limit 
of decimals in the result files. 


time since free-fall 


3 2 1 0 = 
1.5 T T T T T T T 


ro 
oo 


time 


Figure 5. An example of a time symmetric orbit. The orbit is calculated 
first from free-fall position to escape and then correctly back in time to 
the original free-fall situation which is found at t = 3.6 and indicated by 
an arrow. Finally, the integration was continued until the system breaks up 
again at t = 7.2. 


which is based on the equations of motion. However, in most cases 
in our numerical calculations this is not the case even though the 
free-fall moment can be recovered: most orbits actually forget their 
past after the free-fall moment producing an asymmetrical time evo- 
lution. This difference is illustrated with examples in Fig. 5 for time 
symmetric and in Fig. 6 for time asymmetric orbits. We define a 
new quantity £ as time ratio of the recovered part to the orbit’s time 
of disruption fio. If E£ = 1 the orbit is exactly the same in the past 
and in future from the free-fall moment. In Figs 7 and 8 the value 
of & is presented for each orbit identified by its initial position on 
the homology map. The orbits marked by pink colour represent the 
totally symmetrical evolution while other colours indicate binned 
values of £ for asymmetric orbits. In red pixels the length of the 
symmetric part is between 0.45 and 0.60 of tot. In the yellow pixels 
symmetry covers 0.6 and 0.9 of fot and at the blue pixels the orbit 
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Figure 6. An example of an orbit which is asymmetric with respect to time. 
This orbit is also traced back to its free-fall moment at t = 15.8 (marked 
with long arrow) but the motion after this differs from the previous history at 
t = 18.6 (marked with two arrows) leading in this case dramatically different 
breakup. 
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Figure 7. The homology map for most recovered orbits. The structure of 
the map is very similar to the one obtained by ejection times in Fig. 1. Here 
the colour indicates £, the amount of time symmetry, for each orbit. Notice 
how this value tend to increase when initial configuration depart from the 
resonance bands. 


is over 90 but less than 100 per cent symmetric. The white areas — 
distributed quite equally everywhere on the map outside the reso- 
nances — loose their symmetry faster than within 0.45 of tot or they 
represent cases where the orbit has not been recovered. The bins for 
& in Fig. 7 are based on the distribution of orbits in Fig. 9 which we 
will discuss later. 

According to Fig. 7 one will find the exact time symmetric be- 
haviour within resonance areas where orbits have relatively short 
lifetimes. However, is there any relation between the lifetime of 
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Figure 8. The homology map for recovered orbits with £ < 0.45. These 
orbits do not form any grouping as seen for higher values of £ in Fig. 7. 
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Figure 9. The symmetricity, £, of the orbit as a function of the lifetime. The 
functional forms £4 (ftot) and €2(ttor) for two distinctive populations of time 
asymmetric orbits are also shown (see text). Totally time symmetric orbits 
form short horizontal line at £ = 1. 


the asymmetric orbits and their sensitivity to the numerical errors? 
This question is answered in Fig. 9 where we present the amount 
of symmetry, £, as a function of disruption time of the system. It 
appears that the population of asymmetric orbits is divided into two, 
partially overlapping, groups. Both of these groups have distinctive 
functional form £; (tot) = aito + bi and these sets appear to have 
quite sharp edges on the minimum side. Population I (£ = 0.45) 
seems to be direct continuum of time symmetric orbits ( = 1) 
while the lower edge of the curve has approximately parameters 
a, = 0.5, = 0.7 and bı = 0.47. These orbits appear to be located 
just outside the resonance areas on the homology map and the value 
of £ increases when the initial configuration moves farther from the 
resonance bands. 

Most orbits, however, belong to Population II (0 < € < 1) which 
is quite equally distributed over the homology map, clearly further 
from the resonance bands. The lower limit of £ in this population 
obeys approximate relation a) = 2.5,n. = 0.9 and b) = 0. On 
the homology map (Fig. 8) this population is strongly intermixed 
with non-recovered orbits but the latter ones do not show any clear 
functional form for E = £ (t,t). 

We have found out that even in quite a simple dynamical system 
such the classical equal-mass free-fall three-body problem, actually 
the arrow of time shows up. Considering the analytical point of 
view the autonomous form (i.e. no explicit time dependency) of the 
equations of motion and the lack of dissipation suggest that this 
arrow would not be there if we could solve the problem exactly. 
However, in nature totally unperturbed systems do not exist and 
even a minute perturbation, just like a small computing error, is 
adequate to define a definite sense of time. 

How do we find out what is the forward sense of time in the three- 
body problem? The answer is provided by the concept of entropy 
again. Not the classical Boltzmann entropy, but Kolmogorov-Sinai 
entropy (Kolmogorov 1958; Sinai 1959). This concept of entropy 
applies to small systems. It is related to the change of the volume 
of phase space occupied by systems which are initially very close 
to each other. It measures the loss of information about the initial 
conditions or the growth of disorder as time goes on. We have 
previously shown that Kolmogorov-Sinai entropy increases with 
time in the three-body problem, and that it thus tells us what is 
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future, and what is past (Heinämäki et al. 1999). The Kolmogorov— 
Sinai entropy increases towards future time. 

There are other examples of determining the sense of time in the 
three-body problem. In three-body scattering a single body collides 
with a binary, and after some orbital evolution, one of the bodies 
escapes and leaves again a binary. Let us say that the initial binary 
has zero eccentricity. The final binary has an eccentricity e that 
follows the distribution f(e) = 2e (Jeans 1919; Heggie 1975). Thus 
the probability that the final binary escapes with zero eccentricity 
is zero. In a diagram depicting such a scattering event (see e.g. Hut 
& Bahcall 1983) it is immediately clear which is the direction of 
orbital motion even if it is not indicated. Although some of our orbits 
perform time symmetric past and future they resemble a minority in 
our ensemble and it is somewhat reasonable to consider the arrow 
of time in the universe as a statistical quantity. One may also ask if 
these time symmetric orbits really appear at all in true nature where 
small perturbations are always present. 

The microscopic world consists of many small subsystems 
similar to the classical three-body system. It is our conjecture 
that our conclusions would apply also to such systems, and that 
Kolmogorov—Sinai entropy would actually provide the arrow of 
time for the microworld. 

There have previously been a number of papers considering a re- 
lation between the Kolmogorov-—Sinai entropy and the Boltzmann 
thermodynamic entropy. These studies which have been done in 
other fields of physics. They have either tried to address the rela- 
tion between these two quantities, albeit such relation still remains 
unknown (see e.g. Latora & Baranger 1999; Baranger, Latora & 
Rapisarda 2001), or attempted to find a direct equivalence be- 
tween the dynamical instability and the thermodynamic irreversibil- 
ity (Elskens & Prigogine 1986; Yukalov 2003; Ruiz & Tsallis 2007). 
If these relations indeed exist, they endorse our conclusion. 


4 NOISE ESTIMATES 


Since noise estimates are difficult to perform in these kinds of 
systems. We followed the simpler practice of reducing the integra- 
tion step by a factor of 2. We noticed only a slight sharpening in 
the boundaries, but no qualitative differences nor any unexpected 
quantitative differences. 

A second strong test we have performed is for the robustness 
of the method. We ran a naive simulation with non-regularized 
equations of motion but with a short variable step of integration. 
The code is very different in nature from our main code. Yet the 
results are very similar showing the main resonant structures. On 
small scales similar fractal-like patterns appear as with our main 
code. 


5 DISCUSSION 


The results from Anosova (1986) and Heinämäki et al. (1999) sug- 
gested that a homology triangle-like representation is an ideal means 
to describe the evolution or the final solutions of a triple system, and 
that it would show structure both due to global and local sensitivity 
to initial conditions. The geometry of escape times which we have 
calculated shows rich structure featuring resonance bands at various 
levels and areas of apparent higher stability as well as areas with 
clearly a fractal structure. 


The time-scales in question for solar mass objects and an initial 
length scale of R = 1 pc correspond to t ~ 1.6 x 10’ yr. The im- 
plication is that triple stellar systems forming initially from distant 
stars may still be in a multiple system after a few tens of million 
years (the blue areas in Fig. 1 imply time-scales >1.6 x 108 yr). 
On the other hand, one would expect stars being ejected by this 
method from young systems to form loose associations at the latest 
in about a 10 Myr. If the stars were closer (R = 0.1 pc) then the 
breakup times would be a factor of 30 shorter. On the other hand, if 
the stars were not initially at rest then these time-scales are likely to 
be somewhat longer based on previous work (see e.g. Saslaw et al. 
1974; Valtonen & Karttunen 2006). 

Galaxies at distances of ~1 Mpc, and with masses M ~ 50 x 
10° Mo have a typical time-scale of 10!°5 yr, so these are expected 
to be still these initial orbits, even if they were near the resonance 
configuration. 


ACKNOWLEDGMENTS 


This work has been funded by Finnish Cultural Foundation and 
Jenny and Antti Wihuri Foundation. 


REFERENCES 


Aarseth S. J., Anosova J. P., Orlov V. V., Szebehely V. G., 1994, Celest. 
Mech., 58, 1 

Agekian T. A., Anosova J. P., 1967, SvA, 44, 1261 

Anosova J. P., 1986, Ap&SS, 124, 217 

Baranger M., Latora V., Rapisarda A., 2001, Chaos Solutions Fractals, 13, 
471 

Elskens Y., Prigogine I., 1986, Proc. Natl. Acad. Sci. USA, 83, 5756 

Hausdorff F., 1919, Math. Ann., 79, 157 

Heggie D. C., 1975, MNRAS, 173, 729 

Heinämäki P., Lehto H., Valtonen M., Chernin A. D., 1999, MNRAS, 310, 
811 

Huntingford G. W. B., 1961, J. R. Anthropol. Inst. GB Irel., 91, 251 

Hut P., Bahcall J., 1983, ApJ, 268, 319 

Jeans J., 1919, MNRAS, 79, 408 

Kolmogorov A. N., 1958, Dokl. Akad. Nauk SSSR, 119, 861 

Latora V., Baranger M., 1999, Phys. Rev. Lett., 82, 520 

Mikkola S., 1994, MNRAS, 269, 127 

Mikkola S., Aarseth S. J., 1990, Celest. Mech., 47, 375 

Mikkola S., Aarseth S. J., 1993, Celest. Mech., 57, 439 

Mikkola S., Aarseth S. J., 1996, Celest. Mech., 64, 197 

Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745 

Poincaré H., 1892, Les Méthodes Nouvelles de la Méchanique Celeste. 
Gauthier- Villars, Paris 

Prigogine I., Stengers I., 1984, Order Out of Chaos. Bantam Books, New 
York 

Ruiz G., Tsallis C., 2007, Physica A, 386, 720 

Saslaw W. C., Valtonen M. J., Aarseth S. J., 1974, ApJ, 190, 253 

Sinai Ya. G., 1959, Dokl. Akad. Nauk SSSR, 125, 1200 

Umehara H., Tanikawa K., 2000, Celest. Mech., 76, 187 

Valtonen M. J., Karttunen H., 2006, The Three-Body Problem. Cambridge 
Univ. Press, Cambridge 

Yukalov V. I., 2003, Physica A, 320, 149 


This paper has been typeset from a TEX/IATEX file prepared by the author. 


© 2008 The Authors. Journal compilation © 2008 RAS, MNRAS 388, 965-970 


egzoz Jequiaides y, uo IsenB Aq Z7ZEEG6/S96/E/88E/elolWe/SesuW/WOod'dno‘oslwepedce//:sdyjy wo} popeojumoq 


